unit NumerousMethods;

{
  Численные методы
  (c) Keeper
  (r) Василькова Н. Н.
}

{
  PasCalc
  Version 1.2 2011-05-26

  Copyright (C) 2011 KEEPER <yar.keeper@gmail.com>

  This file is part of the PasCalc Calculator

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.

  See the file LICENSE.TXT, included in this distribution,
  for details about the copyright.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
}

interface

uses
  Math, MatrixOperations;

type

  //Таблично заданная функция вида
  //|x0|x1|..|xn|
  //|y0|y1|..|yn|
  TFunction = array [0..1] of TVector;  //0 - args; 1 - values

  { Обобщённый интерполятор. Содержит абстрактные методы }
  TInterpolator = class
  private
    fSourcePoints: TFunction; //Исходная таблица значений функции
    fInterpolant: TVector; //Вектор с результатом интерполяции
    fUnknownPoints: TVector; //Вектор с неизвестными точками
    procedure SetSourcePoints(const AValue: TFunction);
    procedure SetUnknownPoints(const AValue: TVector);
    //Подразумевается, что большая часть необходимых вычислений
    //выполняется после задания таблицы значений функции. Данный метод
    //ВСЕГДА должен быть переопределён!
    procedure AfterSetSourcePoints; virtual; abstract;
  public
    //Конструктору передаётся таблица значений функции SourceTable
    constructor Create(SourceTable: TFunction);
    property SourcePoints: TFunction read fSourcePoints write SetSourcePoints;
    property UnknownPoints: TVector read fUnknownPoints write SetUnknownPoints;
    property InterpResult: TVector read fInterpolant;
    //Заполняет вектор fInterpolant рассчитанными по fUnknownPoints значениями
    procedure FindUnknownPoints;
    //Метод поиска значения в конкретной точке Х. Уникальна для конкретного
    //интерполятора. Всегда должен быть переопределен!
    function GetPoint(const X: extended): Extended; virtual; abstract;
  end;

  { Интерполятор методом Лагранжа }
  TLagrangeInterpolator = class(TInterpolator)
  private
    fQDenominators: TVector;
    fQVector: TVector;
    procedure FormQDenominatorVector;
    function GetQDenominator(const idx: integer): extended;
    function GetQNumenator(const idx: integer; const X: extended): extended;
    procedure AfterSetSourcePoints; override;
    function GetQPoint(i: integer): extended;
  public
    property GetQ[i: integer]: extended read GetQPoint;
    function GetPoint(const X: extended): Extended; override;
    function QCount: integer;
  end;

  TNewtonInterpolator = class(TInterpolator)
  private
    fYDeltas: TVector;
    fHValue: extended;
    fXDeltas: TVector;
    Matrix: TMatrix;
    procedure CalcH;
    procedure CalcXDeltas(const X: extended);
    procedure CreateYDeltas;
    procedure AfterSetSourcePoints; override;
  public
    property YDeltas: TMatrix read Matrix;
    property H: extended read fHValue;
    function GetPoint(const X: extended): Extended; override;
  end;

  TApproximator = class(TInterpolator)
  private
    fDegree: integer;
    fDegreesVector: TVector;
    EquationSystem: TMatrix;
    fSolutionVector: TVector;
    fXYSum: TVector;
    procedure FormXYSumVector;
    procedure FormEquationSystem;
    procedure AfterSetSourcePoints; override;
    procedure GetDegreesVector;
    procedure SetDegree(const Value: integer);
    function GetPolynomCoefficient(i: integer): extended;
  public
    constructor Create(SourceTable: TFunction; PolynomDegree: integer);
    property Matr: TMatrix read EquationSystem;
    property Degrees: TVector read fDegreesVector;
    property XYsum: TVector read fXYSum;
    property PolynomCoefficient[i: integer]: extended read GetPolynomCoefficient;
    function CoefficientsCount: integer;
    property Degree: integer read fDegree write SetDegree;
    function GetPoint(const X: extended): Extended; override;
  end;

implementation

{ TLagrangeInterpolator }

procedure TLagrangeInterpolator.AfterSetSourcePoints;
begin
  FormQDenominatorVector;
end;

procedure TLagrangeInterpolator.FormQDenominatorVector;
var
  i: integer;
begin
  SetLength(fQDenominators, Length(fSourcePoints[0]));
  for i := 0 to high(fQDenominators) do
    fQDenominators[i] := GetQDenominator(i);
end;

constructor TInterpolator.Create(SourceTable: TFunction);
begin
  SetSourcePoints(SourceTable);
end;

procedure TInterpolator.FindUnknownPoints;
var
  i: integer;
begin
  SetLength(fInterpolant, Length(fUnknownPoints));
  for I := 0 to High(fUnknownPoints) do
    fInterpolant[i] := GetPoint(fUnknownPoints[i]);
end;

procedure TInterpolator.SetSourcePoints(const AValue: TFunction);
begin
  fSourcePoints := AValue;
  AfterSetSourcePoints;
end;

procedure TInterpolator.SetUnknownPoints(const AValue: TVector);
begin
  fUnknownPoints := AValue;
end;

function TLagrangeInterpolator.GetQDenominator(const idx: integer): extended;
var
  i: integer;
begin
  result := 1;
  for i := 0 to idx - 1 do
    result := result * (fSourcePoints[0, idx] - fSourcePoints[0, i]);
  for i := idx + 1 to high(fSourcePoints[0]) do
    result := result * (fSourcePoints[0, idx] - fSourcePoints[0, i]);
end;

function TLagrangeInterpolator.GetQNumenator(const idx: integer;
  const X: extended): extended;
var
  i: integer;
begin
  result := 1;
  for i := 0 to idx - 1 do
    result := result * (X - fSourcePoints[0, i]);
  for i := idx + 1 to high(fSourcePoints[0]) do
    result := result * (X - fSourcePoints[0, i]);
end;

function TLagrangeInterpolator.GetPoint(const X: extended): Extended;
var
  i: integer;
begin
  result := 0;
  SetLength(fQVector, Length(fQDenominators));
  for i := 0 to high(fSourcePoints[0]) do
  begin
    fQVector[i] := (GetQNumenator(i, X))/(fQDenominators[i]);
    result := result + (fSourcePoints[1, i] * (fQVector[i]));
  end;
end;

function TLagrangeInterpolator.QCount: integer;
begin
  result := Length(fQVector);
end;

function TLagrangeInterpolator.GetQPoint(i: integer): extended;
begin
  result := fQVector[i];
end;

{ TNewtonInterpolator }

procedure TNewtonInterpolator.AfterSetSourcePoints;
begin
  CreateYDeltas;
  CalcH;
end;

procedure TNewtonInterpolator.CalcH;
begin
  fHValue := fSourcePoints[0, 1] - fSourcePoints[0, 0];
end;

procedure TNewtonInterpolator.CalcXDeltas(const X: extended);
var
  i: integer;
  Val: extended;
begin
  SetLength(fXDeltas, High(fSourcePoints[0]));
  val := 1;
  for I := 0 to High(fXDeltas) do
  begin
    val := val * (X - fSourcePoints[0, i]);
    fXdeltas[i] := val;
  end;
end;

procedure TNewtonInterpolator.CreateYDeltas;
var
  i, j: integer;
begin
  SetLength(Matrix, Length(fSourcePoints[0]), Length(fSourcePoints[0]));
  for I := 0 to high(Matrix) do
    Matrix[i, 0] := fSourcePoints[1, i];
  for i := 1 to High(Matrix) do
    for j := 0 to High(Matrix) - i do
      Matrix[j, i] := Matrix[j+1, i-1] - Matrix[j, i-1];
  SetLength(fYDeltas, High(Matrix));
  for I := 1 to Length(fYDeltas) do
    fYDeltas[i - 1] := Matrix[0, i];
end;

function TNewtonInterpolator.GetPoint(const X: extended): Extended;

function factorial(const n: integer): extended;
var
  i: integer;
begin
  result := 1;
  for I := 2 to n do
    result := result * i;
end;

var
  I: Integer;
  Numenator, Denominator: extended;
begin
  CalcXDeltas(X);
  result := fSourcePoints[1, 0];
  for I := 1 to High(fSourcePoints[0]) do
  begin
    Numenator := fYDeltas[i-1]*fXDeltas[i-1];
    Denominator := factorial(i)*power(fHValue, i);
    result := result + Numenator/Denominator;
  end;
end;

{ TApproximator }

procedure TApproximator.AfterSetSourcePoints;
begin
  FormEquationSystem;
  fSolutionVector := GausseMethod(EquationSystem);
end;

function TApproximator.CoefficientsCount: integer;
begin
  result := Length(fSolutionVector);
end;

constructor TApproximator.Create(SourceTable: TFunction; PolynomDegree: integer);
begin
  fDegree := PolynomDegree;
  inherited Create(SourceTable);
end;

procedure TApproximator.FormEquationSystem;
var
  i, j: integer;
begin
  SetLength(EquationSystem, Degree + 1, Degree + 2);
  GetDegreesVector;
  for i := 0 to high(EquationSystem) do
    for j := 0 to high(EquationSystem[i]) - 1 do
      EquationSystem[i, j] := fDegreesVector[j + i];
  FormXYSumVector;
  for i := 0 to high(EquationSystem) do
    EquationSystem[i, high(EquationSystem[0])] := fXYSum[i];
end;

procedure TApproximator.FormXYSumVector;
var
  i, j: integer;
  sum: extended;
begin
  SetLength(fXYSum, Length(EquationSystem));
  for i := 0 to high(fXYSum) do
  begin
    sum := 0;
    for j := 0 to high(fSourcePoints[0]) do
      sum := sum + (fSourcePoints[1, j] * Power(fSourcePoints[0, j], i));
    fXYSum[i] := sum;
  end;
end;

procedure TApproximator.GetDegreesVector;
var
  i, j: integer;
  sum: extended;
begin
  SetLength(fDegreesVector, 2*Degree + 1);
  for i := 0 to high(fDegreesVector) do
  begin
    sum := 0;
    for j := 0 to high(FSourcePoints[0]) do
      sum := sum + power(fSourcePoints[0, j], i);
      fDegreesVector[i] := sum;
  end;
end;

function TApproximator.GetPoint(const X: extended): Extended;
var
  i: integer;
begin
  result := 0;
  for i := 0 to high(fSolutionVector) do
    result := result + fSolutionVector[i] * power(X, i);
end;

function TApproximator.GetPolynomCoefficient(i: integer): extended;
begin
  result := fSolutionVector[i];
end;

procedure TApproximator.SetDegree(const Value: integer);
begin
  fDegree := Value;
end;

end.

